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Abstract We Consider models of thin disks (with and without bulges) in the 
Bekenstein-Milgrom formulation of MOND as a modification of Newtonian 
gravity. Analytic solutions are found for the full gravitational fields of Kuzmin 
disks, and of disk-plus-bulge generalizations of them. For all these models a 
CTn ■ simple relation between the MOND potential field, i/j, and the Newtonian 

potential, (p N , holds everywhere outside the disk: ^(|V , i/>|/a )V'i/ ; = Vtp N . We 
give exact expressions for the rotation curves for these models. We also find 
that this algebraic relation is a very good approximation for exponential disks. 
The algebraic relation outside the disk is then extended into the disk to derive 
an improved approximation for the MOND rotation curve of disk galaxies that 



requires only knowledge of the Newtonian curve and the surface density. 

I. Introduction 

There are two extreme interpretations of the modified Newtonian dynam- 
ics (MOND). One of these views MOND as a modification of inertia (Milgrom 
1983a, 1994a): Gravitational fields of massive bodies remain Newtonian, but 
the equation of motion of a particle in the field is superseded by a MOND 
equation of motion. In this paper, however, we concentrate on the Bekenstein- 
Milgrom (BM) formulation of MOND (Bekenstein and Milgrom 1984, here- 
after BM), which is an embodiment of MOND as a modification of gravity, 
leaving the Newtonian law of motion intact. The standard Poisson equation 
for the Newtonian gravitational potential, (p N , (V ■ Vf N = AnxGp) induced by 
a mass density p(R) is replaced by 

V-[M(|W|/a )VV>]=47rGp, (1) 

with a the acceleration constant of MOND. This non-linear equation is hardly 



amenable to analytic solution beyond the simple cases of configurations with 
one-dimensional symmetry. 

It would be very useful, for example, to have exact, or even approximate, 
analytic solutions for the gravitational field of model disk galaxies on which 
various ideas can be tested. Some problems whose study may benefit from the 
availability of such solutions are, for example, that of polar rings, and that of 
the motion and fate (disruption, capture etc.) of dwarf companions moving 
in the field of a mother galaxy. 

Even more central is the problem of calculating the MOND rotation 
curves of disk galaxies. In formulations of MOND based on modification of 
inertia, the velocity on a circular orbit of radius r in the plane of disk galaxies 
is given exactly by 

fj,(a/a )a = a N , (2) 

where a = v 2 /r, and a N is the Newtonian acceleration at r (Milgrom 1994a). 
This has been the standard expression for calculating MOND rotation curves 
(e.g. Kent 1987, Milgrom 1988, Begeman, Broeils, and sanders 1991). It is 
not exact in the Bekenstein-Milgrom formulation, and had had the status of 
only an approximation before the work of Milgrom 1994a. 

Here we describe a class of disk-galaxy models for which exact solutions of 
the MOND field equation are presented; this is done in §111. We also find (see 
§IV) that an approximate analytic solution applies for a wider class of models, 
and we suggest a way to predict the adequacy of such an approximation, by 
studying only the Newtonian solution for the mass distribution (§11). In §V 
we describe an approximation for the rotation curve in the BM formulation- 
based, like relation (2) , only on knowledge of the Newtonian acceleration, 
but which is, generally, a better approximation. In §VI we mention further 
possible developments. 

II. An algebraic relation between the Newtonian and MOND 
fields 

Subtracting the Poisson equation from the MOND equation (1) we get 
V • [M(|VV>|/a )VV> - V<p N ] = 0, (3) 

by which the expression in parentheses is some curl field. For configurations 
with one-dimensional symmetry (spherical, cylindrical, or plane) the curl field 
must vanish, and thus the MOND field is related to the Newtonian field by 
the algebraic relation 

M|V^|/a )W> = V<^. (4) 
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This affords a simple solution of the MOND problem by solving first the Pois- 
son equation for ip N , then inverting eq. (4) to get the MOND field. Relation 
(4) does not follow from the MOND equation (1) , but the inverse is correct 
as the latter is just the divergence of the former. 

We begin by asking whether such a relation may hold for more general 
mass distributions, at least approximately. Because the function \i that ap- 
pears in MOND is such that I(x) = xfj,(x) is monotonic, and varies between 
and oo as x does so, I(x) is invertible on the positive real axis. Equation (4) 
is thus equivalent to 

Vt/> = v(\V<p N \/a )V<p N , (5) 

where u(y) = I~ 1 (y)/y. A potential tj; that satisfies this equation exists if and 
only if the curl of the right-hand side vanishes, or, in other terms 

V|V</? W | x Vtp N =0 (6) 

(as v' 7^ 0). This, in turn, is tantamount to |V</?jv| being some function of ip N 

|V^|=/(^)- (7) 

We find then a necessary and sufficient condition for eq. (4) to hold for some 
t/j and some \x (with // ^ 0) ; the condition is expressed solely in terms of the 
Newtonian field of the given mass distribution. By eq. (4) the equipotentials 
for i/j and tp N coincide, and ip is thus a function of cp N . 

A potential ip that satisfies eq. (4) in some domain D is the MOND 
solution to the problem only if if; also satisfies the correct boundary conditions. 
If the sphere at infinity is part of the boundary of D, then i/j automatically 
satisfies the correct boundary condition there. The same is true of the jump 
condition across a thin sheet of mass. If ip N satisfies it than a if) that obeys 
eq. (4) (outside the mass sheet) satisfies the correct jump condition as well. 

Concentrate now on mass distributions that model disk galaxies: An ax- 
isymmetric distribution, symmetric also about a mid-plane, made of a thin 
disk of surface density S(r), and some bulge-like component. By the above 
arguments if ip satisfies eq. (4) everywhere outside the disk it is the MOND 
solution of the problem: the boundary conditions are now satisfied automat- 
ically by a solution of eq. (4) . At infinity V?/> — »■ (MGa ) l/2 R/R 2 , and just 
outside the surface of the disk 

M|W>|/a )cW; = ±27r£(r), (8) 

where d n is te normal component of the gradient. 

To assess the applicability of the algebraic relation for a given configura- 
tion we only have to find the Newtonian potential, and plot |V<^jv| vs - ^jv for 
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points outside the disk. If the points fall on a line, i.e. if |Vv?at| is a function 
of ip N (a highly non-generic case) than Vift as given by eq. (4) is the exact 
MOND acceleration field outside the disk. If |Vy?jv| an d <Pn are correlated 
with only a little scattering, equation (4) gives a good approximation to the 
MOND field (see §IV for examples). 

III. Exact solutions for Kuzmin disks and generalizations thereof 

The two-parameter family of Kuzmin disks is described by a Newtonian 
gravitational potential 

^ K = -MG/[r 2 + (\z\ + h) 2 ] 1/2 (9) 

(see e.g. Binney and Tremaine 1987), where we use cylindrical coordinates 
r, z. The potential above the disk (z > 0) is that of a point mass M placed on 
the lower z-axis at —h = (0,0, —h); the potential below the disk is produced 
by the same mass oppositely placed at h. The surface density, T, K (r), matches 
the jump in the ^-gradient of the potential. 



E*(r) = (2nG)-^ K 



Oz 



Mh/2Ti(r 2 + h 2 ) 3/2 . (10) 



2=cr 



Everywhere outside the disk the equipotential surfaces are concentric spheres 
centered at ±h. Equations (6) (7) are thus satisfied (in this case |V<^at| = 
(p 2 N /MG), and, by the arguments of §11, the exact MOND solution for Kuzmin 
disks is given, outside the disk, by the algebraic relation eq. (4) . Thus, outside 
the disk, 

<f= -Vip-a I~ 1 (g N /a )g N /g N , (11) 

where 

g N = -MG(R±h)/\R±h\ 3 (12) 

is the Newtonian acceleration field above (+), and below (— ) the disk. The 
MOND solution, above the disk, is simply that of a point mass located at —h. 
For very-low-acceleration Kuzmin disks (with MG/h 2 <C a ) we have 
n{x) ps x, so I^ 1 {x) Ri x l/2 . Then, the MOND potential is 

ij K » (MGa ) 1/2 ln[r 2 + (\z\ + h) 2 } 1 ' 2 , (13) 

which can be obtaind by direct integration of eq. (11) . The MOND rotation 
curve of a Kuzmin disk is, by eq. (11) 

v 2 (r) =a o I- 1 ig N (r,0+)/a o ]r 2 /(r 2 + h 2 y/ 2 , (14) 
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where <7jv(r, + ) = MG/(r 2 + h 2 ). It is clear then that we can write 

v 2 (r)=v 2 ca r ] (Cu), (15) 

where v^ = (MGa ) 1/4 is the asymptotic rotational speed, ( = MG/h 2 a is 
a measure of how deep in the MOND regime we are, and u = r/h. If we take, 
for instance, /j(x) = x/(l + x 2 ) l/2 , then I-\y) = [y 2 /2 + (y 2 + y^/A) 1 ' 2 } 1 ' 2 , 
and 



v (r) = Voo 



v 2 



1 + u 



2 



c 2 
1 + 



4(1 + u 



2^2 



1/2 



1/2 



+ WT^> ■ (16) 



In the limit of very-low-acceleration disks, ( — > 0, one has, independently 
of the exact form of /J,(x), 

v 2 (r)=vlu 2 /(l + u 2 ), (17) 

as in this limit I~ 1 (y) = y 1//2 . 

Milgrom (1994b) has proved a virial-like relation for self-gravitating, low- 
acceleration systems, in the BM formulation. For thin disks this relation reads 
(Milgrom 1994b) 

-M 3/2 {Ga ) 1/2 = / 2nrZ(r)v 2 (r) dr, (18) 
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where v(r) is the circular rotation speed; it can be readily verified to hold for 
the pair of £(r) and v(r) given by eqs. (10) and (17) respectively. 

Kuzmin disks may be generalized into a family of disk-plus-bulge mod- 
els that are exactly solvable in MOND. These may be generated in several 
equivalent ways. 

For example, beginning with the Newtonian potential of a Kuzmin disk 
ip K (R) we define a new mass distribution whose Newtonian potential is 

V = U{ VK ). (19) 

We choose U such that U(x) — > x for x — > 0; thus, at spatial infinity <p has 
the same behavior as ipx, and satisfies the correct boundary behavior for a 
potential of a mass M. The potential cp is produced, outside the disk, by a 
mass distribution 

p{R) = (47rG)- 1 W = {AkG)- 1 U"{vk){Vvk) 2 , (20) 

where we have made use of the fact that V 2 (px = 0. From eq. (20) , the 
equidensity surfaces coincide with the equipotential surfaces (common to <p 
and <Pk), because (Vip K ) 2 is a function of <p K . 



In addition, a disk is needed at z = 0, with surface density 



E(r) = (2nG) (h 



z = 



+ 



U'[<p K (r,0)]E K (r), (21) 



with (fxif, 0) = —MG/(r 2 + h 2 ) 1//2 . For p to be non-negative we must have 
U > 0; thus, U is an increasing function. Since the maximum value of (p K is 
0, and there U = 1, we have U < 1 everywhere, or S(r) < T, K (r). The total 
mass (bulge plus disk) contained within an equipotential surface (p K is 

M{ip K ) = U'{ip K )M K {ip K ), (22) 

where M K is the mass within tp K for the generating Kuzmin disk. This can 
be seen by applying the Gauss theorem to the equipotential surface. 

All the potentials <p defined by eq. (19) satisfy eqs. (6) (7) because <p K 
does. Thus the algebraic relation (4) gives the MOND solutions for all these 
model galaxies in term of the Newtonian field V</? = U (<Pk)V(Pk- 

A different approach, which generates the same family of solvable mod- 
els, starts with some spherical density distribution that is centered at —h: 
p(R) = /3(g), q = [r 2 + (\z\ + h) 2 } 1 ' 2 . Take the MOND potential in the z > 
region to coincide with that of p(R). In the z < region the potential is 
defined symmetrically. For spherical systems the MOND field is related to 
the Newtonian field by the algebraic relation (4) . Thus, this is also the 
case for the model under construction. The "bulge" density that produces 
the potential is just the part of the spherical density distribution p(R) that 
is above the mid-plane; we can dictate it at will. A disk with surface den- 
sity S(r) must supplement the bulge to match the jump in the z-gradient. If 
M(q) = J 9 47rA 2 p(A)<iA is the spherical mass within distance q = (r 2 + h 2 ) 1/2 
from the centre of p(R), then 

S(r) is just the surface density of a Kuzmin disk with the same h and a 
mass equal to the total spherical mass within the sphere going through the 
point at r on the disk. Comparing with eq. (21) we find the corresponding 
U'((p K ) = M(g)/M(oo), with g = -MG/tp K . 

A third approach, which we shall not detail here, is to start with the 
MOND potential for the Kuzmin disk, tfj K , and construct new potentials ip = 

S(iPk). 

We reiterate that in all the above models, the bulge equidensity surfaces 
coincide with equipotential surfaces of the model. This means that we can 
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readily construct, for the bulge, distribution functions with isotropic velocity 
distributions. These are of the form f(E), with E = v 2 /2 + ifi(r), for which 
p{f) = Jd 3 vf(E) =F[iP(r)]. 

IV. Some other disk-galaxy models 

How good an approximation is the algebraic relation in general? There 
clearly are disk models for which it fails rankly. Consider, for example a disk 
whose surface density vanishes at the centre. Then, |Vy?jv| vanishes both near 
the centre, and at infinity, while the potential, which vanishes at infinity, is 
non-zero at the centre. Thus |V<^jv| is anything but a function of <p N , and the 
algebraic approximation must break appreciably. 

We have found that for the very pertinent case of a disk with an expo- 
nential surface-density law, £(r) = T< exp(—r/h), the algebraic approximation 
holds very well. A disk for which it holds less well is the so-called Kalnajs 
disk (characterized by a constant angular velocity on circular orbits inside the 
material disk), whose surface density is E(r) = E [l — (r/h) 2 ] 1 ^ 2 . We now 
discuss these two examples in more detail. 

As explained in §11 , to be able to foretell the quality of the algebraic 
approximation for a given disk, it is enough to look at the tightness of the 
relation |V<^jv| vs. ifN- In Fig- 1 we show this relation for the above two 
surface-density distributions, as obtained from numerical calculations using 
a multigrid scheme. The code is capable of solving the nonlinear MOND 
equation, and is described in detail in Brada 1994 (in preparation). For ref- 
erence we also show in Fig. la the numerical results for the Kuzmin disk, 
which show that the numerical scattering about the expected exact relation, 
|V<^at| = <£ 2 N /MG (marked by crosses), is quite negligible (the slight departure 
from the exact relation is numerical, and stem from the cutoff in the disk at 
the end of the mesh) . The correlation for the exponential disk (Fig. lb) is also 
rather tight (but does not follow the asymptotic relation). We thus expect the 
algebraic approximation to be rather good for the MOND field, for all values 
of the mean acceleration. We plot in Fig. 2 the relative departure, 5, from the 
algebraic relation: 



\V<p 



(24) 



AM 



For a very-low-acceleration Kuzmin disk we see that 5 ~ everywhere, as 
expected. For an exponential disk in the same limit (E <C a /G), we see that 
\S\ <C 1 everywhere outside the disk, in keeping with the tight |Vy?Ar| vs. <p N 
relation. For the Kalnajs disk, we see in Fig. lc that the |V<^jv| vs. tp N 
relation has rather more scattering, and indeed the plot of 5, shown in Fig. 2c 
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(again for E <C a /G), evinces more substantial departure from the algebraic 
approximation. An exponential disk with a hole within one-and-a-half scale 
lengths is even a more extreme case described in Figs. Id and 2d. 

V. Rotation curves based on the algebraic approximation 

If the algebraic relation (4) holds outside the disk it cannot be correct 
in the mid-plane of the thin disk; so, we cannot use eq. (2) [n(a/a )a = a N ] 
to obtain the rotation curve of the model galaxy. Rather, we have to follow 
the following procedure: We need the radial acceleration, a r , in the mid-plane 
of the disk. As the acceleration component parallel to the disk is continuous 
across the thin disk, a r is the same as a+, the radial acceleration just outside 
the disk. This can be obtained from the algebraic relation in terms of the total 
Newtonian acceleration just outside the disk a%, and its radial component. 
The latter can again be equated to its value in the mid-plane of the disk (as 
it too is continuous), and so we obtain 

ix{a + /a ) /i[I- 1 (a%/a )]' 

To complete the expression we express a^ in terms of Newtonian radial ac- 
celeration in the mid-plane, directly related to the Newtonian rotation curve 
a N = WIn + (27rGE) 2 ] 1/2 . The MOND rotation curve is thus given by a sim- 
ple function the corresponding Newtonian quantity. The correction to eq. (2) 
involves the addition of the 2-kGT, term in the argument of J -1 in eq. (25) . 

Relation (2) was found numerically (Milgrom 1986) to constitute a good 
approximation for a large class of bulge-plus-disk galaxy models, but, as we 
said, it is not exact in the BM formulation, even for configurations for which 
it is correct outside the disk. For example, for the low-surface-density Kuzmin 
disk, relation (2) gives for the rotation speed 

v 2 (r) = (MGa ) 1/2 r 3/2 /(r 2 + h 2 f'\ (26) 

to be compared with the somewhat different exact expression (17) [ where 
the r dependence is r 2 /(r 2 + h 2 )]. This latter is obtained from eq. (25) . 

We suggest that Eq. (25) is, generically, a better approximation for the 
rotation curve of disk galaxies in the BM formulation than is eq. (2) even 
when the algebraic approximation is not so good outside the disk (see some 
examples below); it is as easy to apply as the latter. [We remind the reader 
that in the formulation of MOND as a modification of inertia (Milgrom 1994a) 
relation (2) gives the rotation curve exactly] 

We give in Fig. 3 three rotation curves for each of a few galaxy models. 
The galaxy models presented are the bare Kuzmin disk, a bare exponential 



disk, and a Kalnajs disk, all in the deep MOND limit. We give the exact rota- 
tion curve calculated numerically, the curve calculated from the approximate 
expression (2) , and that calculated from what we propose as an improved 
approximation (25) . We expect the performance of eq. (2) to be the worst 
for pure disks in the deep MOND limit; adding a spherical component, and/or 
going nearer the Newtonian regime can only improve its performance (but not 
that of approximation (25) ). 

Interestingly, the MOND rotation curve for a Kalnajs disk seems to be 
given exactly by v oc r-as in the Newtonian limit-no matter how deep in the 
MOND regime we are. We do not yet understand the origin of this behavior. 
Once this is taken as fact, the proportionality factor, i.e., the constant angular 
velocity, O, may be calculated, for very-low-acceleration disks from the virial 
relation (18) for disks to get O 2 = 5 ■ 3 _3/2 (27rS G r a o ) 1/2 //i. compared with 
the Newtonian angular velocity which is Q 2 N = n 2 GT, /2h. 

VI. Discussion 

We have described models of disk galaxies for which exact solutions of 
the Bekenstein-Milgrom field equation can be obtained in the form of a simple 
algebraic relation between the MOND solution, and the Newtonian field of the 
same mass distribution. This relation holds approximately for a wider class 
of configurations, which include exponential disks. We have given a simple 
criterion to assess the validity, or near validity of this relation; the use of 
this criterion assumes knowledge of the Newtonian field, (p N , only: it requires 
that |V</9jv| be tightly correlated with <p N outside the disk. We have also 
suggested an improved approximation-inspired by the above approximation- 
for calculating rotation curves in the BM formulation. 

When accuracy beyond the algebraic approximation is needed it may serve 
a first approximation around which we can linearize the MOND equation in 
the small increment. We may, for instance, proceed as follows: Suppose that 
the |V<^jv| vs - <Pn has some scattering but we can reasonably define a mean 
relation |V<^iv| ~ /(^jv)- The acceleration field that is derived from the 
algebraic relation, which is to serve as our zeroth order approximation, is not 
derivable from a potential, in general. So, it is more convenient to work with 
accelerations, not with potentials. Define then 

<f = M (|VV>|/a )VV>, (27) 

which is inverted, as in eq. (5) , to give 

VV> = v(q/a )q, (28) 
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where q = \q\. The algebraic relation would equate q to Vtp N , but we now 
write 

q = V(p N + ff, (29) 

with ff = |V(^at|5 a curl field which is assumed to be small compared with 
S/ifN, and which we shall treat to first order. By the MOND equation we now 
have 

V-rf=0, (30) 

and from eq (28) 

V x [v(q/a )q] = 0. (31) 

Equations (30) (31) are equivalent to the original MOND equation for the 
potential ip (see Milgrom 1986). We now substitute eq. (29) in eq. (31) and 
take the first order in ff (noting that V|V</?jv| x V^w, which measures the 
departure from the algebraic relation, is also first order) to get 

V x ff + ve x [V(e • ff) - f'{^ N )ff\ = z>V|Vv?tv| x e, (32) 

Here e(r) = —Vtp N /\Vtp N \ is a unit vector in the direction of the local New- 
tonian acceleration, and v{f) is the logarithmic derivative of v calculated at 
\Vip N \/a (z> = 1 in the deep MOND limit). The linear equations (30) (32) 
determine ff. 

Our construction of the solvable disk models began with a known MOND 
solution which does not involve a disk, such as a point mass, or, in general, a 
spherical mass distribution. We than place that mass distribution anywhere 
relative to the z = plane; then we take the MOND potential, above the plane 
only, to be that of the mass in question, defining the potential below the plane 
as the mirror image of the one above. The disk is then found which matches 
the jump of the ^-gradient of the potential; hence a family of solvable disk 
models is born. Clearly, we may start with any axisymmetric mass distribution 
for which the MOND solution is known analytically or numerically-not just 
a spherical one-and get a new family of disk models. Such initial non-disk 
MOND solutions can be found by starting from a potential, then calculating 
the density distribution from eq. (1) , making sure that the resulting p is 
positive everywhere, and is otherwise reasonable. 
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Figure captions 

Figure 1. Plots of | V<p N \ vs. tp N for a Kuzmin (a), exponential (b), and Kalnajs 
(c) disks, and for an exponential disk cutoff below one-and-a-half scale lengths. 
The crosses mark the relation |V<^jv| — (Pn/MG. 

Figure 2. A plot of <5-a measure of the departure from the algebraic relation- 
for the four disks as in Fig. 1. 

Figure 3. The rotation curves for the first three disk models of Fig. 1: The 
line is the exact curve; triangles and squares mark the curves calculated, re- 
spectively, by the approximations (2) and (29) . 



12 



